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Abstract. We report large-scale simulations of the three-dimensional Edwards- Anderson 
Ising spin glass system using the recently introduced multi-overlap Monte Carlo algorithm. In 
this approach the temperature is fixed and two replica are coupled through a weight factor 
such that a broad distribution of the Parisi overlap parameter q is achieved. Canonical 
expectation values for the entire q-range (multi-overlap) follow by reweighting. We present 
an analysis of the performance of the algorithm and in particular discuss results on spin glass 
free-energy barriers which are hard to obtain with conventional algorithms. In addition we 
discuss the non-trivial scaling behavior of the canonical (/-distributions in the broken phase. 
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1 Introduction 

The intuitive picture for spin glasses and other systems with conflicting constraints, for 
reviews see |l[ , is that there exists a large number of degenerate thermodynamic states 
with the same macroscopic properties but with different microscopic configurations. 
These states are separated by free-energy barriers in phase space, caused by disorder 
and frustration. Experimentally this is supported by long characteristic times found 
in phenomena like the measurement of the remanent magnetization in typical spin 
glasses such as, e.g., (Feo.isNio.ss^sPieBeAla g]. However, one difficulty of the theory 
of spin glasses is to give a precise meaning to this classification: No explicit order 
parameter exists which allows to exhibit the barriers. The way out of this problem 
appears to use an implicit parametrization, the Parisi overlap parameter q, which 
allows to visualize at least some of them Q. Calculations of the thus encountered 
barriers in q are of major interest. For instance, it is unclear whether the degenerate 
thermodynamic states are separated by infinite barriers or whether this is just an 
artifact of mean-field theory. 

Before performing numerical calculations of these barriers, one of the questions which 
ought to be addressed is "What are suitable weight factors for the problem?" The 
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weight factor of canonical Monte Carlo (MC) simulations is exp(—/3E), where E is 
the energy of the configuration to be updated and (3 is the inverse temperature in nat- 
ural units. The Metropolis and other MC methods generate canonical configurations 
through a Markov process. However, by their very definition the rare-event states 
associated with the free-energy barriers are suppressed in such an ensemble. Now, it 
became widely recognized in recent years that MC simulations with a-priori unknown 
weight factors, like for instance the inverse spectral density l/n(E), are also feasible 
and deserve to be considered, for reviews see j3). Along such lines progress has been 
made by exploring [4-7] innovative weighting methods for the spin glass problem. 
The main idea of the studies [4-7] is to avoid getting stuck in metastable low-energy 
states by using a Markov process which samples the ordered as well as the disordered 
regions of configuration space in one run. Refreshing the system in the disordered 
phase clearly benefits the simulations, but the performance has remained below early 
expectation. One reason appears to be that the direct (i.e. ignoring the dynamics of 
the system) barrier weights are not affected, such that the simulation slows down due 
to the tree-like structure of the low-energy spin glass states, see Ref. || for a detailed 
discussion. 

In Ref. |j] two of the authors introduced a method which focuses directly on enhanc- 
ing the probability for sampling the barrier regions in the Parisi overlap parameter. 
Relying on this method, we have meanwhile performed large-scale simulations of the 
three-dimensional Edwards- Anderson Ising spin glass. First results from this novel 
investigation are reported in the following. The next section briefly introduces the 
spin glass model and is followed by an introduction to the method of || in section 3. 
Technical details of the computer simulation are described in section 4 and some of 
the obtained physical results are summarized in section 5, followed by conclusions in 
the final section 6. 

2 Model 

We focus on the three-dimensional (3c?) Edwards- Anderson Ising (EAI) spin glass on 
a simple cubic lattice of size N — L 3 with periodic boundary conditions. It is widely 
considered to be the simplest model to exhibit realistic spin glass behavior and has 
been the testing ground of Refs. [4-7]. The energy is given by 



where the sum is over nearest-neighbour sites and the Ising spins s, and Sk take 
values ±1. The exchange coupling constants Jik arc quenched random variables which 
are chosen to be ±1 with equal probability. Each fixed assignment of the exchange 
coupling constants Jik defines a realization of the system, and all physical results 
refer to an average over many such realizations. Early MC simulations of the 3d EAI 
model located the freezing temperature at (3 C rj 0.9. For a concise review, see Ref. [Q. 
Recent, very high-statistics canonical simulations |l0j estimate [3 C — 0.901 ± 0.034, 
and were interpreted [jjj to improve the evidence in support of a second-order phase 
transition at j3 c . 

The Parisi overlap parameter is defined as Q 
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where the spins = ±1 and sf = ±1 correspond to two independent copies (replica) 
of the same realization (defined by its couplings Jik), each with its own time evolution 
in the MC simulation (realized by different random numbers). 

3 Multi-overlap algorithm 

The method of Ref. || simulates two replica of the same realization in one computer 
run. In another context this has before been done in (H). Our basic observation, 
closely related to multicanonical methods || || , is that one does still control canonical 
expectation values at temperature when one simulates with a weight function 



Of particular interest is to determine S(q) recursively such that the g-distribution 
becomes uniform in q ("multi-overlap"), and the interpretation of S(q) being the 
microcanonical entropy of the Parisi order parameter. Hence, although an explicit 
order parameter does not exist, an approach very similar to the multimagnetical [ fj"l| 
(which is an highly efficient way to sample interface barriers for ferromagnets) exists 
herewith. 

As a measure of the performance of the algorithm we recorded the average number 
of sweeps that are necessary to perform tunneling events of the form 

(q = 0) — ► (q = ±1) and back. 

In the following we refer to this quantity as autocorrelation time. 

4 Simulation 

The multi-overlap algorithm is thus particularly designed for simulations of the inter- 
esting region below the freezing temperature where Pi(q), the canonical probability 
density of q for realization i (additional dependence on lattice size and temperature 
is implicit), exhibits a rather complicated behavior. The shape of Pi{q) ranges from 
a simple double-peak structure to involved structures of several minima and maxima. 
A few examples encountered in the simulations of the 12 3 system are shown in Fig. |IJ 
This is the situation which is notoriously difficult to study with standard algorithms. 
In our EAI study we therefore focussed on simulations at j3 = 1 > (3 C rj 0.9. We 
investigated lattices of size N = L 3 with L = 4, 6, 8, and 12, using T3E parallel 
computers in Grenoble, Berlin, and Jiilich. We simulated 4096 different realizations 
for the smaller systems up to L = 8, and 512 for L = 12, using the Marsaglia random 
number generator RANMAR for drawing the Jik- As a check we also simulated 
a second set of 4096 realizations for L = 4, 6, and 8, where the Jik were generated 
with a more elaborate version of RANMAR: RANLUX (ll|] with luxury level 4. In 
the simulations themselves we always employed the RANMAR generator due to CPU 
time considerations. 

For each realization the simulation consisted of three steps: 

1. Construction of the weight function (Q). Here we employed an improved variant 
of the accumulative iteration scheme discussed in Ref. IQ. Details will be 
published elsewhere. The iteration was stopped after at least 10 tunnclings 
(between 4 — 20 for L = 12) occurred. 
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Fig. 1: Zoo of probability densities Pi(q) for the 12 3 system at f3 = 1. The subscript 
i = l,..., 512 labels the different realizations. 



2. Equilibration run. This run was to equilibrate the system for given fixed weight 
factors. 

3. Production run. Each production run of data taking was concluded after at least 
20 tunneling events were recorded. To allow for easy (standard) reweighting in 
temperature we stored besides the Parisi overlap parameter also the energies and 
magnetizations of the two replica in a time-series file. Since the corresponding 
autocorrelation times r vary significantly from realization to realization we used 
an adaptive data compression routine to make sure that the spacing between 
measurements was always proportional to r. Using this adaptive scheme we 
recorded for each realization 65536 measurements. 

Let us conclude this section with a few remarks on the performance of the new al- 
gorithm. Fitting the estimates of the mean autocorrelation time [r] av to the form 
ln([T] a v) = a + z ln(iV) gives z = 2.42 ± 0.03. As will be discussed below, the implied 
improvement with respect to barrier calculations is huge. Nevertheless the slowing 
down is quite off from the theoretical optimum, which is z = 1 for multicanonical sim- 
ulations j|. One reason seems to be that we want the g-distribution to be uniform 
in the whole admissible interval q S [— 1, +1], including the region of |g| ss 1 that is 
strongly correlated with ground states and hence difficult to reach by local updates, 
see for instance ||. Being content with a smaller region (like the region between the 
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Fig. 2: Integrated probability density F of canonical tunneling barrier heights at /3 = 1. 



two outmost maxima of Pi(q)) is expected to give further improvements of the tun- 
neling performance. By also monitoring the encountered minimum, maximum and 
median tunneling times we observed that the mean values are systematically larger 
than the median, what means that the tunneling distribution has a rather long tail 
towards large tunneling times. On the other hand, the effect is not severely hindering 
our multi-overlap simulations: For the lattice sizes L = 4 to 8 the worst behaved real- 
ization took never more than 1% of the entire computer time. Indicating a remarkable 
increase of complexity, this amount was about 12% for L = 12. 



5 Results 

The data created in this way allows us to calculate a number of physically interesting 
quantities. Let us first concentrate on the canonical potential barriers in q which were 
in Ref. || defined as 

-Aq 

B>= H max [1, Pi(q)/Pi(q + Aq)} , (4) 

<?=-! 

where Aq is the step-size in q. For the double-peak situations of first-order phase 
transitions Jll| Eq. (|) simplifies to B t = pm^/p™* w here i^ max is the absolute 
maximum and P™ 111 is the absolute minimum (for ferromagnets at q = 0) of the 
probability density Pi{q). Our definition generalizes to the situation where several 
minima and maxima occur due to disorder and frustration. When evaluating (jij) from 
numerical data for Pi(q) some care is needed to avoid contributions from statistical 
fluctuations of Pi (q) . 

Graphically, our values for the Bi are presented in Fig. Q It comes as a surprise 
that the finite-size dependence of the distributions is very weak. To study this issue 
further, we have compiled in Table [j] for each lattice size the following informations 
about our potential barrier results: largest and second largest values B max and i?2, 
median values P mc d and mean values [P] av with their statistical error bars. From this 
table it becomes obvious, why this investigation could not be performed using canon- 
ical methods to which in this context also multicanonical simulations and enlarged 
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Table 1: Canonical potential barriers: maximum (and its contribution to the mean 
in %), second largest value, the median and its jackknife error bar, the mean and its 
(standard) error bar. 



L 


^max 


B 2 






4 


8.76E07 (62%) 


6.01E06 


11.05 ±0.24 


(1.72 ± 1.08)E04 


G 


3.73E08 (41%) 


3.63E08 


12.83 ±0.45 


(1.12±0.64)E05 


8 


1.54E09 (77%) 


1.41E08 


14.52 ±0.26 


(2.44 ± 1.89)E05 


12 


9.14E07 (73%) 


1.78E06 


13.29 ±0.88 


(1.53 ± 1.44)E05 



ensembles belong, as their weights for those barriers are still canonical. For these 
methods the slowing down is proportional to the average barrier height [-B] av , which 
is already large for L = 4, about 17 thousand, and increases to about 250 thousand 
for L — 8. The subsequent decrease to about 150 thousand for L = 12 has to be 
attributed to the smaller number of realizations in this case: Comparison of the ex- 
trem values makes only sense when the numbers of realizations match. On L — 4 and 
L = 6 systems we have performed a number of (very long) canonical simulations to 
estimate the proportionality constant between barrier height and improvement due to 
our multi-g simulations. Using these result, we estimate that with our computer pro- 
gram a canonical MC calculation of the worst L = 12 barrier alone, the one reported 
in Table [j], would take about 1000 years on a 500 MHz Alpha processor. 
The reader may be puzzled by the very large error bars assigned to the mean barrier 
values. Their explanation is: The entire mean value is dominated by the largest 
barrier, which contributes between 41% (L — 6) and 77% (L — 8), see B meix in the 
second column of Table |l|. Besides B max , the second largest value B 2 is listed in 
the third column. It may be remarked that most of these worst case barriers exhibit 
simple double-peak behavior. One lesson from these numbers is that very few of the 
realizations are responsible for the collapse of canonical simulation methods. 
Typical realizations, described by the median results of Table [l], have much smaller 
tunneling barriers. They turn out to be quite insensitive to the lattice size. This 
result of an almost constant typical tunneling barrier is consistent with the fact that 
our tunneling times are rather far apart from their theoretical optimum: Still other 
reasons than overlap barriers have to be responsible. Therefore, one may question the 
apparently accepted opinion that these are the typical barriers which are primarily 
responsible for the severe slowing down of canonical MC simulations. 
Even though less detailed than the barrier data also the averaged canonical probability 
densities P(q) = [Pi{q)]&v contain important information on the system. Due to the 
average over all realizations the P(q) exhibit a double-peak shape; see Figs. ||^. In 
principle the data stored in the time-series files allow canonical reweighting in a p- 
range around the simulation point. For disordered systems, however, where many 
realizations have to be reweighted, great care is needed when estimating the valid 
reweighting range. 

By analyzing the spin glass susceptibility, xsg = -^[(g 2 )] av , we obtained in Ref. (s) 
the best finite-size scaling fit xsg °c IP'" at /? = 0.88 with 7/1/ = 2.37(4) and a 
goodness-of-fit parameter Q = 0.25. This was corroborated by the curves of the 
Binder parameter, g = (l/2)(3 — [(q 4 )]a,v /[{q 2 )]1 v ), which merge around (3 = 0.89. In 
the low-temperature phase (f3 > f3 c ) the curves for different lattice sizes seem to fall 
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Fig. 3: Finite-size scaling plot for P(q) at the transition temperature. 



on top of each other, but our error bars were still too large to draw a firm conclusion 
from this quantity. 

Close to the transition temperature T c one expects that the probability densities P(q) 
for different lattice sizes satisfy the finite-size scaling relation 

P(q)=L^P(L^q,L 1 / v (T-T c )) , (5) 

where P is a scaling function, and (3 and v are the critical exponents of the order 
parameter and correlation length, respectively. They are related to 7/1/ = 2 — 77 by 
the standard scaling relation f3/v — (d — j/i/)/2. Using our estimate for 7/1/ w 2.37 
we thus obtain the value (3/v ~ 0.317. Figure || shows the probability densities 
P{q) of the large-scale simulations, reweighted to the transition temperature T c and 
rescaled according to (||), using the above exponent estimates. We see that the data 
for L = 4, 6, and 8 fall almost perfectly onto a common master curve, while the 
L = 12 data obviously show some deviations, in particular close to the peak. The 
overall appearance of Fig. ||, therefore, is actually worse than that of the corresponding 
plot in Ref. || which is based on a much smaller exploratory data set. We strongly 
suspect that the reason for the deviations of the L — 12 curve can be traced back to 
problems with the admissible reweighting range for our largest system size. Since for 
quenched, disordered systems the reweighting range for averaged quantities depends 
on the common overlap of the reweighting ranges for the individual realizations, it is 
indeed conceivable that this problem does show up more dramatically for the larger 
sample of realizations. Given the very good data collapse for the smaller lattices, and 
with this explanation in mind, we feel that also our new results are compatible with 
the findings of Ref.JlOfl. We are currently trying to improve the test of the finite-size 
scaling prediction (p) by redoing the simulations closer to j3 c . Narrowing the g-range 
may allow to simulate lattices of size L = 16 and beyond. 

The observation that the Binder parameter curves do not splay out at low temper- 
atures suggests that, similar to the 2d XY model, the correlation length £ may be 
infinite - or at least larger than the simulated lattice sizes. Generalizing the second 
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argument of F to L/( and assuming L/£ = 0, one thus expects that the P(q) should 
scale also below T c . With our simulation set-up, the most reliable test of this conjec- 
ture can be done at the simulation point T = 1 » 0.88T C , since then no temperature 
reweighting is involved. The raw data for P{q) at T = 1 are shown in Fig. Notice 
that -P(O) is slightly growing with increasing system size. By adjusting the only free 
parameter, j3/v — 0.255, we obtain the finite-size scaling plot in Fig. [3| which shows 
a very clear data collapse onto a single master curve for all lattice sizes. Moreover, 
if we reweight our data to the even lower temperatures T = 1/1.1 (= 0.80T C ) and 
T = 1/1.2 (= 0.73T C ) we still find a reasonable data collapse, see Fig. [| but here 
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Fig. 6: Finite-size scaling plot for P(q) at T « 0.8T C . 



again extreme care is needed to control the reliable reweighting range. 
Of course, since our lattice sizes are relatively small, we cannot conclude that the 
correlation length is infinite below T c . If £ is large but finite, it is conceivable that 
we observe an effective scaling behavior as long as £ > L. We may conclude, that the 
correlation length is unusually large (£ > 12) down to 0.73T C . 



6 Conclusions 

For the 3d EAI model at f3 = 1 we have performed a high-statistics calculation for 
probability distributions depending on the Parisi overlap parameter q. The results 
for free-energy barriers in q became feasible by using g-dependent (multi-overlap) 
weight factors in large-scale MC simulations. Although the tunneling performance is 
not optimal, the method opens new horizons for spin glass simulations. Using slight 
modifications of the method (like narrowing the g-range, including a magnetic field, 
etc.) will allow us to extend our investigation into various interesting directions, like 
an improved study of the thermodynamic limit at and below the freezing point, a 
study of the 4d EAI spin glass model, or e-physics, where one studies the influence of 
an interaction term pl| eJ^iLi s \ s 1 — q in the Hamiltonian ([[]). In multi-overlap 
simulations we can obtain expectation values for arbitrary e-values by reweighting. 
Physically most interesting is the case where a non-zero magnetic field is combined 
with a non-zero e- value. 
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